Accepted by the Astrophysical Journal 



FIRST YEAR WILKINSON MICROWAVE ANISOTROPY PROBE (WMAP) 
OBSERVATIONS: TESTS OF GAUSSIANITY 

E. Komatsu ^, A. Kogut ^, M. R. Nolta ^, C. L. Bennett ^, M. Halpern ^, G. Hinshaw ^, N. 
Jarosik ^ M. Limon S. S. Meyer ^ L. Page ^ D. N. Spergel ^, G. S. Tucker 3-6,8^ l_ Verde 2.9, 

E. Wollack 3, E. L. Wright ^° 

komatsu@astro.princeton.edu 
ABSTRACT 

We present Hmits to the amphtude of non-Gaussian primordial fluctuations in the 
WMAP 1-year cosmic microwave background sky maps. A non-linear coupling param- 
eter, /nl, characterizes the amplitude of a quadratic term in the primordial potential. 
We use two statistics: one is a cubic statistic which measures phase correlations of 
temperature fluctuations after combining all configurations of the angular bispectrum. 
The other uses the Minkowski functionals to measure the morphology of the sky maps. 
Both methods find the WMAP data consistent with Gaussian primordial fluctuations 
and establish limits, —58 < /nl < 134, at 95% confidence. There is no significant 
frequency or scale dependence of /nl- The WMAP limit is 30 times better than COBE, 
and validates that the power spectrum can fully characterize statistical properties of 
CMB anisotropy in the WMAP data to high degree of accuracy. Our results also 
validate the use of a Gaussian theory for predicting the abundance of clusters in the 
local universe. We detect a point-source contribution to the bispectrum at 41 GHz, 
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bsrc = (9-5 ± 4.4) X 10~^ /xK^ sr^, which gives a power spectrum from point sources of 

Csrc = (15 =b 6) X 10~^ //K^ sr in thermodynamic temperature units. This value agrees 
well with independent estimates of source number counts and the power spectrum at 
41 GHz, indicating that 6src directly measures residual source contributions. 

Subject headings: cosmic microwave background — cosmology: observations — early 
universe — galaxies: clusters: general — large-scale structure of universe 

1. INTRODUCTION 

The Gaussianity of the primordial fluctuations is a key assumption of modern cosmology, 
motivated by simple models of inflation. Statistical properties of the primordial fluctuations are 
closely related to those of the cosmic microwave background (CMB) radiation anisotropy; thus, a 
measurement of non- Gaussianity of the CMB is a direct test of the inflation paradigm. If CMB 
anisotropy is Gaussian, then the angular power spectrum fully specifles the statistical properties. 
Recently, Acquaviva et al. (2002) and Maldacena (2002) have calculated second-order perturbations 
during inflation to show that simple models based upon a slowly-rolling scalar field cannot generate 
detectable non-Gaussianity. Their conclusions are consistent with previous work (Salopck &: Bond 
1990, 1991; Falk et al. 1993; Gangui et al. 1994). Infiation models that have significant non- 
Gaussianity may have some complexity such as non-Gaussian isocurvature fluctuations (Linde & 
Mukhanov 1997; Peebles 1997; Bucher &; Zhu 1997), a scalar-fleld potential with features (Kofman 
et al. 1991; Wang & Kamionkowski 2000), or "curvatons" (Lyth & Wands 2002; Lyth et al. 2002). 
Detection or nondetection of non-Gaussianity thus sheds light on physics of the early universe. 

Many authors have tested Gaussianity of CMB anisotropy on large angular scales (~ 7°) 
(Kogut et al. 1996; Heavens 1998; Schmalzing &; Gorski 1998; Ferreira et al. 1998; Pando et al. 
1998; Bromley & Tegmark 1999; Banday et al. 2000; Contaldi et al. 2000; Mukherjee et al. 2000; 
Magueijo 2000; Novikov et al. 2000; Sandvik & Magueijo 2001; Barreiro et al. 2000; Phillips & 
Kogut 2001; Komatsu et al. 2002; Komatsu 2001; Kunz et al. 2001; Aghanim et al. 2001; Cayon 
et al. 2002), on intermediate scales (~ 1°) (Park et al. 2001; Shandarin et al. 2002), and on small 
scales (~ 10') (Wu et al. 2001; Santos et al. 2002; Polenta et al. 2002). So far there is no evidence 
for signiflcant cosmological non-Gaussianity. 

Most of the previous work only tested the consistency between the CMB data and simulated 
Gaussian realizations without having physically motivated non-Gaussian models. They did not, 
therefore, consider quantitative constraints on the amplitude of possible non-Gaussian signals al- 
lowed by the data. On the other hand, Komatsu et al. (2002), Santos et al. (2002), and Cayon 
et al. (2002) derived constraints on a parameter characterizing the amplitude of primordial non- 
Gaussianity inspired by inflation models. The former and the latter approaches are conceptually 
different; the former does not address how Gaussian the CMB data are or the physical implica- 
tion of the results. In this paper, we adopt the latter approach, and constrain the amplitude of 
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primordial non-Gaussianity in the WMAP 1-year sky maps. 

Some previous work all had roughly similar sensitivity to non-Gaussian CMB anisotropy at 
different angular scales, because the number of independent pixels in the maps are similar, i.e., 
~ 4000 - 6000 for COBE (Bennett ct al. 1996), QMASK (Xu et al. 2002), and MAXIMA (Hanany 
et al. 2000) sky maps. Polenta et al. (2002) used about 4 x 10^ pixels from the BOOMERanG map 
(de Bernardis et al. 2000), but found no evidence for non-Gaussianity. The WMAP provides about 
2.4 X 10^ pixels (outside the KpO cut) uncontaminated by the Galactic emission (Bennett et al. 
2003a), achieving more than one order of magnitude improvement in sensitivity to non-Gaussian 
CMB anisotropy. 

This paper is organized as follows. In § 2, we describe our methods for measuring the primordial 
non-Gaussianity using the cubic (bispectrum) statistics and the Minkowski functionals, and present 
the results of the measurements of the WMAP 1-year sky maps. Implications of the results for 
inflation models and the high-redshift cluster abundance are then presented. In § 3, we apply the 
bispectrum to individual frequency bands to estimate the point-source contribution to the angular 
power spectrum. The results from the WMAP data are then presented, and also comparison among 
different methods. In § 4, we present summary of our results. In Appendix A, we test our cubic 
statistics for the primordial non-Gaussianity using non-Gaussian CMB sky maps directly simulated 
from primordial fluctuations. In Appendix B, we test our cubic statistic for the point sources using 
simulated point-source maps. In Appendix C, we calculate the CMB angular bispectrum generated 
from features in a scalar-field potential. 



2. LIMITS ON PRIMORDIAL NON-GAUSSIANITY 

2.1. The WMAP 1-year Sky Maps 

We use a noise-weighted sum of the Ql, Q2, VI, V2, Wl, W2, W3, and W4 maps. The maps 
are created in the HEALPix format with nside = 512 (Gorski et al. 1998), having the total number 
of pixels of 12 X nside^ = 3, 145, 728. We do not smooth the maps to a common resolution before 
forming the co-added sum. This preserves the independence of noise in neighboring pixels, at the 
cost of complicating the effective window function for the sky signal. We assess the results by 
comparing the WMAP data to Gaussian simulations processed in identical fashion. Each CMB re- 
alization draws a sample from the ACDM cosmology with the power-law primordial power spectrum 
fit to the WMAP data (Hinshaw et al. 2003b; Spergel et al. 2003). The cosmological parameters are 
in Table 1 of Spergel et al. (2003) (we use the best-fit ^^WMAP data only" parameters). We copy 
the CMB realization and smooth each copy with the WMAP beam window functions of the Ql, Q2, 
VI, V2, Wl, W2, W3, and W4 (Page et al. 2003a). We then add independent noise realizations to 
each simulated map, and co-add weighted by Nobs/(^o, where the effective number of observations 
iVobs varies across the sky. The values of the noise variance per N^hs, Cq' ^® tabulated in Table 1 
of (Bennett et al. 2003b). 
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We use the conservative KpO mask to cut the Galactic plane and known point sources, as 
described in Bennett et al. (2003a), retaining 76.8% of the sky (2,414,705 pixels) for the subsequent 
analysis. In total 700 sources are masked on the 85% of the sky outside the Galactic plane in all 
bands; thus, the number density of masked sources is 65.5 sr~^. The Galactic emission outside the 
mask has visible effects on the angular power spectrum (Hinshaw et al. 2003b). Since the Galactic 
emission is highly non-Gaussian, we need to reduce its contribution to our estimators of primordial 
non-Gaussianity. Without foreground correction, both the bispectrum and the Minkowski func- 
tionals find strong non-Gaussian signals. We thus use the forground template correction described 
in section 6 of Bennett et al. (2003c) to reduce foreground emission to negligible levels in Q, V, and 
W bands. The method is termed as an "alternative fitting method", which uses only the Q, V, and 
W band data. The dust component is separately fitted to each band without assuming spectrum of 
the dust emission (3 parameters). We assume that the free-free emission has spectrum, and 

the synchrotron has z/~^-^ spectrum. The amplitude of each component in Q band is then fitted 
across three bands (2 parameters). 



2.2. Methodology 

2.2.1. Model for Primordial Non-Gaussianity 

We measure the amplitude of non-Gaussianity in primordial fluctuations parametrized by a 
non-linear coupling parameter, /nl (Komatsu &; Spergel 2001). This parameter determines the 
amplitude of a quadratic term added to the Bardeen curvature perturbations $ ($h in Bardeen 
(1980)), as 

$(x) = $l(x) + /nl [$E(x) - (^E(x)>] , (1) 

where <I>l are Gaussian linear perturbations with zero mean. Although the form in equation (1) is 
inspired by simple inflation models, the exact predictions from those inflation models are irrelevant 
to our analysis here because the predicted amplitude of /nl is much smaller than our sensitivity; 
however, this parameterization is useful to find quantitative constraints on the amount of non- 
Gaussianity allowed by the CMB data. Equation (1) is general in that /nl parameterizes the 
leading-order non- linear corrections to We discuss the possible scale-dependence in Appendix C. 

Angular bispectrum analyses found |/nl| < 1500 (68%) from the COBE DMR 53-^90 GHz 

coadded map (Komatsu et al. 2002) and |/nl| < 950 (68%) from the MAXIMA sky map (Santos 
et al. 2002). The skewness measured from the DMR map smoothed with filters, called the Spherical 
Mexican Hat wavelets, found |/nl| < 1100 (68%) (Cayon ct al. 2002), although they neglected the 
integrated Sachs-Wolfe effect in the analysis, and therefore underestimated the cosmic variance of 
/nl- BOOMERanG did not measure /nl in their analysis of non-Gaussianity (Polenta et al. 2002). 
The r.m.s. amplitude of $ is given by ($^)^^^ ~ {^1)^^^ (l + /nl (^l))- Since {^'^^^^ measured 
on the COBE scales through the Sachs-Wolfe effect is ($2)^2 ^ ^{/^T'^f''^ /T ~ 3.3 x 10"^ 
(Bennett et al. 1996), one obtains /nl (*l) < 2-5 x lO'^ from the COBE 68% constraints; thus. 
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we already know that the contribution from the non-hnear term to the r.m.s. amphtude is smaUer 
than 0.25%, and that to the power spectrum is smaher than 0.5%. This amphtude is comparable 
to limits on systematic errors of the WMAP power spectrum (Hinshaw et al. 2003a), and needs to 
be constrained better in order to verify the analysis of the power spectrum. 



2.2.2. Method 1: The Angular Bispectrum 

Our first method for measuring /nl is a "cubic statistic" which combines nearly optimally 
all configurations of the angular bispectrum of the primordial non-Gaussianity (Komatsu et al. 
2003). The bispectrum measures phase correlations of field fluctuations. We compute the spherical 
harmonic coefficients a/^ of temperature fluctuations from 

aim = I d^nM(n)^^y4(n), (2) 

where M(n) is a pixel-weighting function. Here, M(n) is the KpO sky cut where M(n) takes 
in the cut region and 1 otherwise. We filter the measured o;^ in Z-space and transform it back to 
compute two new maps, A{r, n) and B(r, n), given by 

A{r,h) = E ^ ^^ai^YUh), (3) 

1=2 m=-l ' 

^max ^ f \hi 

Bir,n) ^ E E ^^aiMn). (4) 

1=2 m=-l ' 

Here Ci = Cibf + A^, where C/ is the CMB anisotropy, N the noise bias, and bi the beam window 
function describing the combined smoothing effects of the beam (Page et al. 2003a) and the finite 
pixel size. The functions a;(r) and Pi(r) are defined by 

= k^dkgTiik)ji{kr), (5) 

(3i{r) = k^dkP{k)gTi{k)jiikr), (6) 

where r is the comoving distance. These two functions constitute the primordial angular bispectrum 
and correspond to ai{r) = f^lbf^{r) and A(r) = bf{r) in the notation of Komatsu & Spergel 
(2001). We compute the radiation transfer function gTi{k) with a code based upon CMBFAST 
(Seljak & Zaldarriaga 1996) for the best-fit cosmological model of the WMAP 1-year data (Spergel 
et al. 2003). We also use the best-fit primordial power spectrum of ^, P{k). We then compute the 
cubic statistic for the primordial non-Gaussianity, iSprim, by integrating the two filtered maps over 
r as (Komatsu et al. 2003) 

5prim = ^ j Airr^dr j ^^(r, ry)B^{r, n), (7) 
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where the angular average is done on the full sky regardless of sky cut, and ms = (47r)~^ / (fnM^{n) 
is the third-order moment of the pixel-weighting function. When the weight is only from a sky cut, 
as is the case here, wc have = /sky, i-e., is the fraction of the sky covered by observations 
(Komatsu et al. 2002). Komatsu et al. (2003) show that i? is a Wiener-filtered map of the underlying 
primordial fluctuations, The other map A combines the bispectrum configurations that are 
sensitive to non-linearity of the form in equation (1). Thus, >Sprim is optimized for measuring the 
skewness of $ and picking out the quadratic term in equation (1). 



Finally, the non-linear coupling parameter /nl is given by 

-1 



/i 



NL 



h<l2<h ^ ^ ^ 



"Sprini) (8) 



where Sf^i^^ is the primordial bispectrum (Komatsu &; Spergel 2001) multiphed by bij^bi^bi^ and 
computed for /nl = 1 and the best-fit cosmological model. This equation is used to measure /nl 
as a function of the maximum multipole Zmax- The statistic <Sprim takes only Ar3/2 operations to 
compute without loss of sensitivity whereas the full bispectrum analysis takes N^^"^ operations. It 
takes about 4 minutes on 16 processors of an SGI Origin 300 to compute /nl from a sky map at 
the highest resolution level, nside = 512. We measure /nl as a function of Zmax- Since there is little 
CMB signal compared with instrumental noise at I > 512, we shall use Zmax = 512 at most; thus, 
nside = 256 is sufficient, speeding up evaluations of /nl by a factor of 8 as the computational time 
scales as (nside)^. The computation takes only 30 seconds at nside = 256. Note that since we are 
eventually fitting for two parameters, /nl and 6src (see sec. 3), we include covariance between these 
two parameters in the analysis. The covariance is, however, small (see Figure 8 in Appendix A). 

While we use uniform weighting for M(n), we could instead weight by the inverse noise variance 
per pixel, M(n) = A^^-'^(n); however, this weighting scheme is sub-optimal at low I where the CMB 
anisotropy dominates over noise so that the uniform weighing is more appropriate. For measuring 
6src) on the other hand, we shall use a slightly modified version of the weighting, as ftgrc comes 
mainly from small angular scales where instrumental noise dominates. 



2.2.3. Method 2: The Minkowski Functionals 

Topology offers another test for non-Gaussian features in the maps, measuring morphological 
structures of ffuctuation fields. The Minkowski functionals (Minkowski 1903; Gott et al. 1990; 
Schmalzing &; Gorski 1998) describe the properties of regions spatially bounded by a set of contours. 
The contours may be specified in terms of fixed temperature thresholds, v = AT/a, where a is 
the standard deviation of the map, or in terms of the area. Parameterization of contours by 
threshold is computationally simpler, while parameterization by area reduces correlations between 
the Minkowski functionals (Shandarin et al. 2002). We use a joint analysis of the three Minkowski 
functionals (area A{u), contour length C(i^), and genus G{i')) explicitly including their covariance; 
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consequently, we work in the simpler threshold parameterization. 

The Minkowski functionals are additive for disjoint regions on the sky and are invariant under 
coordinate transformation and rotation. We approximate each Minkowski functional using the set 
of equal-area pixels hotter or colder than a set of fixed temperature thresholds. The fractional area 

is thus the number of enclosed pixels, Ny, divided by the total number of pixels on the cut sky, 
iVcut- Here is the area of an individual spot, and A is the total area of the pixels outside the cut. 
The contour length 

cw = ^E^^ (10) 

i 

is the total perimeter length of the enclosed regions Pi, while the genus 

= ^(iVhot-A^coid) (11) 

is the number of hot spots, A^hot) minus the number of cold spots, A^coid- We calibrate finite 
pixelization effects by comparing the Minkowski functionals for the WMAP data to Monte Carlo 
simulations. 

The WMAP data are a superposition of sky signal and instrument noise, each with a different 
morphology. The Minkowski functionals transform monotonically (although not linearly) between 
the limiting cases of a sky signal with no noise and a noise map with no sky signal. Unlike spatial 
analyses such as Fourier decomposition, different regions of the sky cannot be weighted by the 
signal-to-noise ratio, nor does the noise "average down" over many pixels. The choice of map 
pixelization becomes a tradeoff between resolution (favoring smaller pixels) versus signal-to-noise 
ratio (favoring larger pixels). We compute the Minkowski functionals at nside = 16 through 256 
(3072 to 786,432 pixels on the full sky). We use the WMAP KpO sky cut to reject pixels near 
the Galactic plane or contaminated by known sources. The cut sky has 1433 pixels at resolution 
nside = 16 and 666,261 pixels at nside = 256. 

We compute the Minkowski functionals at 15 thresholds from —3.5(7 to +3.5a, and compare 
each functional to the simulations using a goodness-of-fit statistic, 

= ^ [-^WMAP ~ (-^sim)]j^i "^uiu^ ["^WMAP ~ (-^sim)]i^2! (12) 

1/11/2 

where -F'-vvmap ^ Minkowski functional from the WMAP data (the index i denotes a kind of 
functional), {F^i^) is the mean from the Monte Carlo simulations, and is the bin-to-bin 

covariance matrix from the simulations. 
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2.3. Monte Carlo Simulations 

Monte Carlo simulations are used to estimate the statistical significance of the non-Gaussian 
signals. One kind of simulation generates Gaussian random realizations of CMB sky maps for 
the angular power spectrum, window functions, and noise properties of the WMAP 1-year data. 
This simulation quantifies the uncertainty arising from Gaussian fields, or the uncertainty in the 
absence of non-Gaussian fluctuations. The other kind generates non-Gaussian CMB sky maps from 
primordial fluctuations of the form of equation (1) (see Appendix A for our method for simulating 
non-Gaussian maps). This simulation quantifies the uncertainty more accurately and consistently 
in the presence of non-Gaussian fluctuations. 

In principle, one should always use the non-Gaussian simulations to characterize the uncer- 
tainty in /nl; however, the uncertainty estimated from the Gaussian reahzations is good approxi- 
mation to that from the non-Gaussian ones as long as |/nl| < 500. Our non-Gaussian simulations 
verify that the distribution of /nl and 6src around the mean is the same for Gaussian and non- 
Gaussian realizations (see Figure 8 in Appendix A for an example of /nl = 100). The Gaussian 
simulations have the advantage of being much faster than the non-Gaussian ones. The former takes 
only a few seconds to simulate one map whereas the latter takes 3 hours on a single processor of 
an SGI Origin 300. Also, simulating non-Gaussian maps at nside = 512 requires 17 GB of physical 
memory. We therefore use Gaussian simulations to estimate the uncertainty in measured /nl and 

^src- 

2.4. Limits to Primordial Non-Gaussianity 

Figure 1 shows /nl measured from the Q-|-V-|-W coaddcd map using the cubic statistic [Eq (8)], 
as a function of the maximum multipole /max- We find the best estimate of /nl = 38 =b 48 (68%) 
for /max = 265. The distribution of /nl is close to a Gaussian, as suggested by Monte Carlo 
simulations (see Figure 8 in Appendix A). The 95% confidence interval is —58 < /nl < 134. 
There is no significant detection of /nl at any angular scale. The r.m.s. error, estimated from 
500 Gaussian simulations, initially decreases as oc /max; although /nl for /max = 265 has a smaller 
error than that for /max = 512 because the latter is dominated by the instrumental noise. Since 
all the pixels outside the cut region are uniformly weighted, the inhomogeneous noise in the map 
(pixels on the ecliptic equator are noisier than those on the north and south poles) is not accounted 
for. This leads to a noisier estimator than a minimum variance estimator. The constraint on 
/nl for /max = 512 will improve with more appropriate pixel-weighting schemes (Heavens 1998; 
Santos et al. 2002). The simple inverse noise {N~^) weighting makes the constraints much worse 
than the uniform weighting, as it increases errors on large angular scales where the CMB signal 
dominates over the instrumental noise. The uniform weighting is thus closer to optimal. Note 
that for the power spectrum, one can simply use the uniform weighting to measure Q at small / 
and the weighting at large /. For the bispectrum, however, this decomposition is not simple. 
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as the bispectrum Bi^i^i^ measures the mode couphng from li to I2 and Z3 and vice versa. This 
property makes it difficult to use different weighting schemes on different angular scales. The first 
column of Table 1 shows /nl measured in Q, V, and W bands separately. There is no a significant 
band-to-band variation, or a significant detection in any band. 

Figure 2 shows the Minkowski functionals at nside = 128 (147,594 high-latitude pixels, each 28' 
in diameter). The gray band shows the 68% confidence region derived from 1000 Gaussian simula- 
tions. Table 2 shows the values [Eq.(12)]. The data are in excellent agreement with the Gaussian 
simulations at all resolutions. The individual Minkowski functionals are highly correlated with each 
other (e.g., Shandarin et al. (2002)). We account for this using a simultaneous analysis of all three 
Minkowski functionals, replacing the 15-element vectors -P^map u and (-Fsim u) equation (12) (the 
index i denotes each Minkowski functional) with 45-element vectors = [F^ , F'^ , F'^]^ = [Area, 
Contour, Genus],/ and using the covariance of this larger vector as derived from the simulations. 
We compute for values /nl = to 1000, comparing the results from WMAP to similar values 
computed from non-Gaussian realizations. Figure 3 shows the result. We find a best-fit value 
/nl = 22 lb 81 (68%), with 95% confidence upper limit /nl < 139, in agreement with the cubic 
statistic. 



2.5. Implications of the WMAP Limits on /nl 

2.5.1. Inflation 

The limits on /nl are consistent with simple infiation models: models based on a slowly rolling 
scalar field typically give |/nl| ~ 10~^ - 10"^ (Salopek & Bond 1990, 1991; Falk et al. 1993; Gangui 
et al. 1994; Acquaviva et al. 2002; Maldacena 2002), three to four orders of magnitude below our 
limits. Measuring /nl at this level is difficult because of the cosmic variance. There are alternative 
models which allow larger amplitudes of non-Gaussiantiy in the primordial fluctuations, which we 
explore below. 

A large /nl may be produced when the following condition is met. Suppose that <1> is given 
by $ = ex where e is a transfer function that converts x to $ and x = x^^^ +x^'^'> + 0(x^^^) denotes 
a fluctuating field expanded into a series of x^'^^ = fiX^^~^^x^^^ with /i = 1. Then, /nl = e~^/2- 
Inflation predicts the amplitude of x^^^ and the form of fi which eventually depends upon the 
scalar field potential; thus, x^*-* would be of order (.ff/mpianck)* {H is the Hubble parameter during 
inflation) for H < mpianck, and the leading order term is eff/mpianck 10~^e. In this way e 
"suppresses" the amplitude of fluctuations, allowing a larger amplitude for i^/rripianck ~ 10~^e~^. 
What does this mean? li H ^ 10~^mpianck) then e ~ 10^^ and /nl ~ 10^/2. The amplitude of /nl 
is thus large enough to detect for /2 > 0.1. This suppression factor, e, seems necessary for one to 
obtain a large /nl in the context of the slow-roll inflation. The suppression also helps us to avoid 
a "fine-tuning problem" of inflation models, as it allows i^/rripianck to be of order slightly less than 
unity (which one might think natural) rather than forcing it to be of order 10~^. 



-10- 



Table 1: The non- linear coupling parameter, the reduced point-source angular bispectrum, and the 
point-source angular power spectrum (positive definite) by frequency band. The errorbars are 68%. 
The tabulated values are for the KpO mask, while the Kp2 mask gives simila r results. 

/nl "src Cgrc 

[10-5 /xK3 sr2] [10-3 /xK2 sr] 



Q 51 ±61 9.5 ±4.4 15 ± 6 

V 42 ±63 1.1 ±1.6 4.5 ±4 

W 37 ±75 0.28 ±1.3 — 

Q+V+W 38 ±48 0.94 ±0.86 — 



Table 2: for Minkowski Functionals' 



nside 


Pixel Diam 


Minkowski 


WMAP 


f{>WMAP)^ 




(deg) 


Functional 


x' 




256 


0.2 


Genus 


15.9 


0.57 


128 


0.5 


Genus 


10.7 


0.79 


64 


0.9 


Genus 


15.7 


0.44 


32 


1.8 


Genus 


18.7 


0.26 


16 


3.7 


Genus 


16.8 


0.22 


256 


0.2 


Contour 


9.9 


0.93 


128 


0.5 


Contour 


9.9 


0.83 


64 


0.9 


Contour 


14.6 


0.54 


32 


1.8 


Contour 


12.8 


0.58 


16 


3.7 


Contour 


11.9 


0.67 


256 


0.2 


Area 


17.4 


0.50 


128 


0.5 


Area 


10.9 


0.74 


64 


0.9 


Area 


11.9 


0.66 


32 


1.8 


Area 


21.9 


0.12 


16 


3.7 


Area 


15.7 


0.33 



"X^ computed using Gaussian simulations. There are 15 degrees of freedom. 
''Fraction of simulations with greater than the value from the WMAP data. 



- 11 - 



1200 
1000 h 

800 

600 



400 



200 h 


-200 h 
-400 



10 



"1 r 



"1 — |— r 



"1 — I — I I I 



Q+V+W coadded map 



f,,= 38±48 for ^,=265 



J I I I I 



100 



max 



J I I I I I I 



1000 



Fig. 1. — The non-linear coupling parameter /nl as a function of the maximum multipole imax, 
measured from the Q+V+W coadded map using the cubic (bispectrum) estimator [Eq. (8)]. The 
best constraint is obtained from Zmax = 265. The distribution is cumulative, so that the error bars 
at each luiax. are not independent. 
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Fig. 2. — The left panels show the Minkowski functionals for WMAP data (filled circles) at nside = 
128 (28' pixels). The gray band shows the 68% confidence interval for the Gaussian Monte Carlo 
simulations. The right panels show the residuals between the mean of the Gaussian simulations 
and the WMAP data. The WMAP data are in excellent agreement with the Gaussian simulations. 
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Fig. 3.— Limits to /nl from fit of tlic WMAP data to the non-Gaussian models [Eq. (1)]. The 
fit is a joint analysis of the three Minkowski functional at 28' pixel resolution. There are 44 degrees 
of freedom. 
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Curvatons proposed by Lyth &; Wands (2002) provide an example of a supression mechanism. 
A curvaton is a scalar field, a, having mass, rricr, that develops fluctuations, Sa, during inflation 
with its energy density, po- — ^(^)) fii^Y compared to that of the inflaton field that drives inflation. 
After inflation ends, radiation is produced as the inflaton decays, generating entropy perturbations 
between a and radiation, S^^ = Spa-/pa — jSp-y/p^- When H decreases to become comparable to 
nia, oscillations of a at the bottom of V{a) give pa — m^a^. In the limit of "cold inflation" for 
which ^p-J p-^i is nearly zero, one finds <S'o-7 — ^ pa I Per — 26a / a + {5a / a)'^ . As long as a survives after 
the production of S^y, the curvature perturbation $ is generated as $ = ^e5cr7 — e[x^^^ + ^(x*-^-*)^] 
where x^^^ = 6a/a (i.e., /2 = ^). The generation of <1> continues until a decays, and <1> is essentially 
determined by a ratio of p^ to the total energy density, Q^, at the time of the decay. Lyth et al. 
(2002) numerically evolved perturbations to find e ~ at the time of the decay. The smaller the 
curvaton energy density is, the less efficient the S^j to $ conversion becomes (or the more efficient 
the supression becomes). The small rig- thus leads to the large /nli as /nl = e~^/2 — j^a^ (i-^-^ 
/nl is always positive in this model). Assuming the curvaton exists and is entirely responsible for 
the observed CMB anisotropy, our limits on /nl imply > 9 x 10^^ at the time of the curvaton 
decay. (However, the lower limit to O^. does not mean that we need the curvatons. This constraint 
makes sense only when the curvaton exists and is entierly producing the observed fluctuations.) 

Features in an inflaton potential can generate significant non-Gaussian fluctuations (Kofman 
et al. 1991; Wang &; Kamionkowski 2000), and it is expected that measurements of non-Gaussianity 
can place constrains on a class of the feature models. In Appendix C, we calculate the angular 
bispectrum from a sudden step in a potential of the form in equation (C2). This step is motivated 
by a class of supergravity models yielding the steps as a consequence of successive spontaneous 
symmetry-breaking phase transitions of many fields coupled to the inflaton (Adams et al. 1997, 
2001). One step generates two distinct regions in I space where |/nl| is very large: a positive /nl 
is predicted at Z < Zf, while a negative /nl at / > Zf, where If is the projected location of the step. 
Our calculations suggest that the two regions are separated in I by less than a factor of 2, and one 
cannot resolve them without knowing If. The average of many I modes further smeares out the 
signals. The averaged /nl thus nearly cancels out to give only small signals, being hidden in our 
constraints in Figure 1. Peiris et al. (2003) argue that some sharp features in the WMAP angular 
power spectrum producing large values may arise from features in the inflaton potential. If this 
is true, then one may be able to see non-Gaussian signals associated with the features by measuring 
the bispectrum at the scales of the sharp features of the power spectrum. 



2.5.2. Massive cluster abundance at high redshift 

Massive halos, like clusters of galaxies at high redshift, are such rare objects in the universe 
that their abundance is sensitive to the presence of non-Gaussianity in the distribution function of 
primordial density fluctuations. Several authors have pointed out the power of the halo abundance 
as a tool for finding primordial non-Gaussianity (Lucchin &; Matarrese 1988; Robinson & Baker 
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2000; Matarrese et al. 2000; Benson et al. 2002); however, the power of this method is extremely 

sensitive to the accuracy of the mass determinations of halos. It is necessary to go to redshifts 
of z > 1 to obtain tight constraints on primordial non-Gaussianity, as constraints from low and 
intermediate redshifts appear to be weak (Koyama et al. 1999; Robinson et al. 2000) (see also 
Figure 4 and 5). Due to the difficulty of measuring the mass of a high-redshift cluster the current 
constraints are not yet conclusive (Willick 2000). The limited number of clusters observed at 
high redshift also limits the current sensitivity. In this section, we translate our constraints on 
/nl from the WMAP 1-year CMB data into the effects on the massive halos in the high-redshift 
universe, showing the extent to which future cluster surveys would see signatures of non-Gaussian 
fluctuations. 

We adopt the method of Matarrese et al. (2000) to calculate the dark-matter halo mass function 
dn/dM for a given /nl, using the ACDM with the running spectral index model best-fit to the 
WMAP data and the large-scale structure data. This set of parameters is best suited for the 
calculations of the cluster abundance. The parameters are in the rightmost column of Table 8 of 
Spergel et al. (2003). We calculate 

dn _ ^PrnO 

dM~ M 

where pmo = 2.775 x 10^^(J7m/i^) Mq Mpc"^ = 3.7 x 10^° Mq Mpc"^ is the present-day mean mass 
density of the universe, P(M, z) is the probability for halos of mass M to collapse at redshift 
and dP/dM is given by 

dM Jo 27r \dM ^ SdM ' ^ ' 

where the angle 9x is given by 6x = X6c + X^Hs/Q, and 6c{z) is the threshold overdensity of spherical 
collapse (Lacey & Cole 1993; Nakamura Sz Suto 1997). The variance of the mass fluctuations as a 
function of z is given by cr^{M, z) = D^{z)a'^{M, 0), where D{z) is the growth factor of linear density 
fluctuations, a'^{M, 0) = dkk-'^Flj{k)A'^{k), A'^{k) = {2TT^)-'^k^P{k) is the dimensionless power 
spectrum of the Bardeen curvature perturbations, FM{k) = g{k)T{k)W{kRM) a filter function, 
g{k) = |(A;/i7o)^J^jnO ^ conversion factor from <I> to density fluctuations, T{k) the transfer function 
of linear density perturbations, W{x) = 3ji{x)/x the spherical top-hat window smoothing density 
fields, and Rm = [3M/(47r/9mo)]^''^ the spherical top-hat radius enclosing a mass M. The skewness 
/Lt3(M, z) = D^{z)iJ,3{M,0), where 

/X3(M,0)=6/nl / -r^FM{ki)A\ki) -lFM{k2)A\k2) / dfiFM{Jkf + kl + 2kik2fi), 

Jo Kl Jo 1^2 Jo ^ 

(15) 

arises from the primordial non-Gaussianity. We use a Monte Carlo integration routine called vegas 
(Press et al. 1992) to evaluate the triple integral in equation (15). It follows from equation (15) 
that a positive /nl gives a positive fi^, positively skewed density fluctuations. Also this dn/dM 
reduces to the Press-Schechter form (Press & Schechter 1974) in the limit of /nl 0. Although the 
Press-Schechter form predicts signiflcantly fewer massive halos than A?^-body simulations (Jenkins 
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et al. 2001), we assume that a predicted ratio of the non-Gaussian dn/dM to the Gaussian dn/dM 

is still reasonably accurate, as the primordial non-Gaussianity does not affect the dynamics of 
halo formations which causes the difference between the Press-Schechter form of dn/dM and the 
A^-body simulations. 

Figure 4 shows the WMAP constraints on the ratio of non-Gaussian dn/dM to the Gaussian 
one, as a function of M and z. We find that the WMAP constraint on /nl strongly limits the 
amplitude of changes in dn/dM due to the non-Gaussianity. At ^ = 0, dn/dM is changed by no 
more than 20% even for 4 x 10^^ Mq clusters. The number of clusters that would be newly found 
at z = 1 for M < 10^^ Mq should be within of the value predicted from the Gaussian theory. 
At z = 3, however, much larger effects are still allowed: dn/dM can be increased by up to a factor 
of 2.5 for 2 X 10^"^ Mq. 

Predictions for actual cluster surveys are made clearer by computing the source number counts 
as a function of z, 



where V{z) is the comoving volume per steradian, and M^j^ is the limiting mass that a survey can 
reach. In practice would depend on z due to, for example, the redshift dimming of X-ray 

surface brightness; however, a constant M^^ turns out to be a good approximation for surveys of 
the Sunyaev-Zel'dovich (SZ) effect (Carlstrom ct al. 2002). Figure 5 shows the ratio for dN/dz as 
a function of z and M\\^. A source-detection sensitivity of S'lim = 0.5 Jy roughly corresponds to 
Miim = 1.4 X 10^^ Mq (Carlstrom et al. 2002), for which dN/dz should follow the prediction of 
the Gaussian theory out to z ~ 1 to within 10%, but dN/dz at z = 3 can be increased by up to a 
factor of 2. As -Miim increases, the impact on dN/dz rapidly increases. 

The SZ angular power spectrum Cf^ is so sensitive to that we can use Cf^ to measure erg 
(Komatsu & Kitayama 1999). The sensitivity arises largely from massive (M > 10^^ Mq) clusters at 
z ~ 1. Prom this fact one might argue that Cf"^ is also sensitive to the primordial non-Gaussianity. 
We use a method of Komatsu & Scljak (2002) with dn/dM replaced by equation (13) to compute 
Cf^ for the WMAP limits on /nl • We find that Cf^ should follow the prediction from the Gaussian 
theory to within 10% for 100 < I < 10000. This is consistent with Cf^ being primarily sensitive to 
halos at z ~ 1, where the effect on dN/dz is not too strong (see Figure 5). Since Cf^ oc (78(Ob/t)^ 
(Komatsu k, Seljak 2002), erg can be determined from Cf^ to within 2% accuracy at a fixed J^b^ 
using the Gaussian theory. The current theoretical uncertainty in the predictions of Cf^ is a factor 
of 2 in Cf^ (10% in as), still much larger than the effect of the non-Gaussianity. 




(16) 



-17- 




Fig. 4. — The limits to the effect of the primordial non-Gaussianity on the dark-matter halo mass 
function dn/dM as a function of z. The shaded area represents the 95% constraint on the ratio of 
the non-Gaussian dn/dM to the Gaussian one. 
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Fig. 5. — The same as figure 4 but for the dark-matter halo number counts dN/dz as a function of 
the hmiting mass of a survey. 
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3. LIMITS TO RESIDUAL POINT SOURCES 

3.1. Point-source Angular Power Spectrum and Bispectrum 

Radio point sources distributed across the sky generate non-Gaussian signals, giving a positive 
bispectrum, bsrc (Komatsu & Spcrgcl 2001). In addition, the point sources contribute significantly 
to the angular spectrum on small angular scales (Tegniark & Efstathiou 1996), contaminating the 
cosmological angular power spectrum. It is thus important to understand how much of the measured 
angular power spectrum is due to sources. We constrain the source contribution to the angular 
power spectrum, Cgrc, by measuring tgrc- Komatsu Sz Spergel (2001) have shown that WMAP can 
detect 6src even after subtracting all (bright) sources detected in the sky maps. Fortunately, there 
is no degeneracy between /nl and ftsro 

as shown later in Appendix A. 

In this section we measure the amplitude of non-Gaussianity from "residual" point sources 
which are fainter than a certain flux threshold, Sc, and left unmasked in the sky maps. The 
bispectrum fegrc is related to the number of sources brighter than Sc per solid angle A'^(> Sc)- 

fSc i^fj rSc 

hUSc) = j^ dS—[gi,y)Sf = -N{>ScM'^)Scf + 3j^ -^N{>SMu)Sf, (17) 

where g{i') is a conversion factor from Jy sr~^ to fxK which depends on observing frequency v as 

g{i^) = (24.76 Jy fiK-^ sr-i)-i[(sinhx/2)/x2]2, x = hu/kBTo ~ z//(56.78 GHz) for Tq = 2.725 K 
(Mather et al. 1999), and dN/dS is the differential source count per solid angle. The residual point 
sources also contribute to the point-source power spectrum Cgrc as 

csrc(^c) = y^ dS—[giu)Sf = -Ni>ScMu)Scf + 2j^ ^N{> SMu)S]' . (18) 

By combining equation (17) and (18) we find a relation between 6src and Cgrc, 

fSc jq 

C^rciSc) = b,,c{ScMl^)Sc]-' + — 6src(5)b(^)-S]-^ (19) 

We can use this equation combined with the measured bgrc as a function of Sc to directly determine 
function of Sc, without relying on any extrapolations. When the source counts obey a 
power-law like dN/dS oc S^, one finds bsrc{S) oc S'^'^^; thus, brighter sources contribute more to 
the integral in equation (19) than fainter ones as long as /3 > —3, which is the case for fluxes of 
interest. Bennett ct al. (2003a) have found /3 = —2.6 it 0.2 for 5 = 2 — 10 Jy in Q band. Below 1 Jy, 
/? becomes even flatter (Toffolatti et al. 1998), implying that one does not have to go down to the 
very faint end to obtain reasonable estimates of the integral. In practice, we use equation (17) with 
N{> S) of the Toffolatti et al. (1998) model (hereafter T98) at 44 GHz to compute bsrc{S < 0.5 Jy), 
inserting it into the integral to avoid missing faint sources and underestimating the integral. 
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3.2. Measurement of the Point-source Angular Bispectrum 

The reduced point-source angular bispectrum, 6src) is measured by a cubic statistic for point 
sources (Komatsu et al. 2003), 

Sps = m,'J ^D\n), (20) 
where the filtered map D{n) is given by 

^(") = E E -^^^imYimin). (21) 

1=2 m=-l 

This statistic is even quicker (~ 100 times) to compute than 5prim (eq-[7]), as it involves only one 
integral over n and only one filtered map. This statistic also retains the same sensitivity to the 
point-source non-Gaussianity as the full bispectrum analysis. The cubic statistic Sps gives bgrc as 

-1 

«Sps, (22) 

where Bf^i^^,^ is the point-source bispectrum for 6src = 1 (Komatsu Sz Spergel 2001) multiplied by 
bi^bi^bi^. While the uniform pixel- weighting outside the Galactic cut was used for /nl, we use 
here M(n) = [ctcmb + -^(n)] ~^ where (t^mb = (47r)~^ X)i(2^ + is the variance of CMB 

anisotropy and A'^(n) is the variance of noise per pixel which varies across the sky. This weighting 
scheme is nearly optimal for measuring 6src as the signal comes from smaller angular scales where 
noise dominates. The factor of c^mb approximately takes into account the non-zero contribution 
to the variance from CMB anisotropy. This weight reduces uncertainties of bgrc by 17%, 23%, and 
31% in Q, V, and W bands, respectively, compared to the uniform weighting. We use the highest 
resolution level, nside = 512, and integrate equation (22) up to Zmax = 1024. In Appendix B, it 
is shown that this estimator is optimal and unbiased as long as very bright sources, which have 
contributions to Ci too large to ignore, are masked. We cannot include Cgrc in the filter, as it is 
what we are trying to measure using bsrc- 

The filled circles in the left panels of Figure 6 represent bgrc measured in Q (top panel) and V 
(bottom panel) band. We have used source masks for various flux cuts, Sq, defined at 4.85 GHz to 
make these measurements. (The masks are made from the GB6-I-PMN 5 GHz source catalogue.) 
We find that 6src increases as Sc- the brighter sources unmasked, the more non-Gaussianity is 
detected. On the other hand one can make predictions for bgic using equation (17) for a given 
A'^(> S). Comparing the measured values of 6src with the predicted values from A'^(> S) of T98 
(dashed lines) at 44 GHz, one finds that the measured values are smaller than the predicted values 
by a factor of 0.65. The solid lines show the predictions multiplied by 0.65. Both errors in the 
T98 predictions and a non-fiat energy spectrum of sources easily cause this factor. (If sources have 
a non-flat spectrum like S oc u"' where a 0, then Sc at Q or V band is different from that at 
4.85 GHz.) Bennett et al. (2003a) find that the majority of the radio sources detected in Q band 
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have a flat spectrum, a = 0.0 ifc 0.2. Our value for the correction factor matches well the one 
obtained from the WMAP source counts for 2 — 10 Jy in Q band (Bennett et al. 2003a). 

Equation (18) combined with the measured ftgrc is used to estimate the point-source angular 
power spectrum Cgrc- The right panels of Figure 6 show the estimated Cgrc as filled circles. These 
estimates agree well with predictions from equation (18) with A^(> S) of T98 multiplied by a factor 
of 0.65 (solid lines). For = 1 Jy at Q band, Cgrc = (19 ± 5) x 10"'^ /liK^ sr, and matches well the 
value estimated from the WMAP source counts at the same flux threshold (Bennett et al. 2003a) , 
which corresponds to the solid lines in the figure. At V band, Cgrc = (5±4) x 10~^ /xK^ sr. Here, the 
hat denotes that these values do not represent Cgrc for the standard source mask used by Hinshaw 
et al. (2003b) for estimating the cosmological angular power spectrum. Since the standard source 
mask is made of several source catalogues with different selection thresholds, it is difficult to clearly 
identify a mask flux cut. We give the standard mask an "effective" flux cut threshold at 4.85 GHz by 
comparing ftgrc measured from the standard source mask (shaded areas in Figure 6; see the second 
column of Table 1 for actual values) with those from the GB6+PMN masks defined at 4.85 GHz. 
The measurements agree when c± 0.75 Jy in Q band. Using this effective threshold, one expects 
Cgrc for the standard source mask as Cgrc = (15 ±6) x 10~^ /xK^ sr in Q band. This value agrees with 
the excess power seen on small angular scales, (15.5 it 1.7) x 10~^ fiK^ sr (Hinshaw et al. 2003b), as 
well as the value extrapolated from the WMAP source counts in Q band, (15.0 ± 1.4) x 10~^ /iK^ sr 
(Bennett et al. 2003a). In V band, c^rc = (4.5 ± 4) x 10"^ //K^ sr. 

The source number counts, angular power spectrum, and bispectrum measure the first-, second- 
, and third-order moments of dN/dS, respectively. The good agreement among these three different 
estimates of Cgrc indicates the validity of the estimate of the effects of the residual point sources in 
Q band. There is no visible contribution to the angular power spectrum from the sources in V and 
W bands. We conclude that our understanding of the amplitude of the residual point sources is 
satisfactory for the analysis of the angular power spectrum not to be contaminated by the sources. 

4. CONCLUSIONS 

We use cubic (bispectrum) statistics and the Minkowski functionals to measure non-Gaussian 
fluctuations in the WMAP 1-year sky maps. The cubic statistic [Eq.(7)] and the Minkowski func- 
tionals place limits on the non-linear coupling parameter /nl, which characterizes the amplitude of 
a quadratic term in the Bardeen curvature perturbations [Eq. (1)]. It is important to remove the 
best-fit foreground templates from the WMAP maps in order to reduce the non-Gaussian Galactic 
foreground emission. The cubic statistic measures phase correlations of temperature fluctuations to 
find the best estimate of /nl from the foreground-removed, weighted average of Q-|-V-|-W maps as 
/nl = 38±48 (68%) and —58 < /nl < 134 (95%). The Minkowski functions measure morphologi- 
cal structures to find /nl = 22 it 81 (68%) and /nl < 139 (95%), in good agreement with the cubic 
statistic. These two completely different statistics give consistent results, validating the robustness 
of our limits. Our limits are 20—30 times better than the previous ones (Komatsu et al. 2002; 



- 22 - 




0.5 1 1.5 0.5 1 1.5 
S^@4. 85GHz [Jy] 

Fig. 6. — The point-source angular bispectrum bgrc and power spectrum Cgrc- The left panels show 
6src in Q band (top panel) and V band (bottom panel). The shaded areas show measurements 
from the WAIAP sky maps with the standard source cut, while the filled circles show those with 
flux thresholds Sc defined at 4.85 GHz. The dashed lines show predictions from equation (17) with 
A^(> S) modeled by Toffolatti et al. (1998), while the solid lines are those multiplied by 0.65 to 
match the WMAP measurements. The right panels show Cgrc- The filled circles are computed from 
the measured 6src substituted into equation (19). The lines are from equation (18). The error bars 
are not independent, because the distribution is cumulative. 
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Santos et al. 2002; Cayon et al. 2002), and constrain the relative contribution from the non-hnear 
term to the r.m.s. amphtude of $ to be smaller than 2 x 10~^ (95%), much smaller than the limits 
on systematic errors in the WMAP observations. This validates that the angular power spectrum 
can fully characterize statistical properties of the WMAP CMB sky maps. We conclude that the 
WMAP 1-year data do not show evidence for significant primordial non-Gaussianity of the form in 
equation (1). Our limits are consistent with predictions from inflation models based upon a slowly 
rolling scalar field, |/nl| = 10"^ - 10"^ The span of all non-Gaussian models, however, is large, 
and there are models which cannot be parametrized by equation (1) (e.g., Bernardeau & Uzan 
(2002b, a)). Other forms such as multi-field inflation models and topological defects will be tested 
in the future. 

The non-Gaussianity also affects the dark-matter halo mass function dn/dM, since the massive 
halos at high redshift are sensitive to changes in the tail of the distribution function of density 
fluctuations. Our limits show that the number of clusters that would be newly found at 2; = 1 for 
M < 10^^ Mq should be within 1^°^ of the value predicted from the Gaussian theory. At higher 
redshifts, however, much larger effects are still allowed. The number counts dN/dz at 2: = 3 with 
the limiting mass of 3 x 10^"^ A/q can be reduced by a factor of 2, or increased by more than a 
factor of 3. Since the SZ angular power spectrum is primarily sensitive to massive halos at z 1 
where the impact of non-Gaussianity is constrained to be within 10%, a measurement of ag from 
the SZ angular power spectrum is changed by no more than 2%. Our results on dn/dM derived in 
this paper should be taken as the current observational limits to non-Gaussian effects on dn/dM. 
In other words, this is the uncertainty that we currently have in dn/dM when the assumption of 
Gaussian fluctuations is relaxed. 

The limits on /nl will improve as the WMAP satellite acquires more data. Monte Carlo 
simulations show that the 4-year data will achieve 95% limit of 80. This value will further improve 
with a more proper pixel-weighting function that becomes the uniform weighting in the signal- 
dominated regime (large angular scales) and becomes the N^^ weighting in the noise-dominated 
regime (small angular scales). There is little hope of testing the expected levels of /nl = 10""^ — 10~^ 
from simple inflation models, but some non-standard models can be excluded. 

We have detected non-Gaussian signals arising from the residual radio point sources left un- 
masked at Q band, characterized by the reduced point-source angular bispectrum 6src = (9.5 it 
4.4) X 10^^ fiK'^ sr^, which, in turn, gives the point-source angular power spectrum Cgrc = (15 ± 
6) X 10^^ fiK^ sr. This value agrees well with those from the source number counts (Bennett 
et al. 2003a) and the angular power spectrum analysis (Hinshaw et al. 2003b), giving us confidence 
on our understanding of the amplitude of the residual point sources. Since bgrc directly measures 
Cgrc without relying on extrapolations, any CMB experiments which suffer from the point-source 
contamination should use 6src to quantify Cgrc to obtain an improved estimate of the CMB angular 
power spectrum for the cosmological-parameter determinations. 

Hinshaw et al. (2003b) found that the best-fit power spectrum to the WMAP temperature 
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data has a relatively large value, corresponding to a chance probability of 3%. While still 

acceptable fit, there may be missing components in the error propagations over the Fisher matrix. 
Since the Fisher matrix is the four-point function of the temperature fluctuations, those missing 
components (e.g., gravitational lensing effects) may not be apparent in the bispectrum, the three- 
point function. The point-source non-Gaussianity contributes to the Fisher matrix by only a 
negligible amount, as it is dominated by the Gaussian instrumental noise. Non-Gaussianity in 
the instrumental noise due to the 1// striping may have additional contributions to the Fisher 
matrix; however, since the Minkowski functionals, which are sensitive to higher-order moments 
of temperature fluctuations and instrumental noise, do not find significant non-Gaussian signals, 
non-Gaussianity in the instrumental noise is constrained to be very small. 

The WMAP mission is made possible by the support of the Office of Space Sciences at NASA 
Headquarters and by the hard and capable work of scores of scientists, engineers, technicians, 
machinists, data analysts, budget analysts, managers, administrative staff, and reviewers. LV is 
supportd by NASA through Chandra Fellowship PF2-30022 issued by the Chandra X-ray Obser- 
vatory center, which is operated by the Smithsonian Astrophysical Observatory for an on behalf of 
NASA under contract NAS8-39073. 



A. SIMULATING CMB SKY MAPS FROM PRIMORDIAL FLUCTUATIONS 

In this appendix, we describe how to simulate CMB sky maps from generic primordial fluc- 
tuations. As a specific example, we choose to use the primordial Bardeen curvature perturbations 
$(x), which generate CMB anisotropy at a given position of the sky Ar(n) as (Komatsu et al. 



where ^im{r) is the harmonic transform of 'l>(x) at a given comoving distance r = |x|, ^imif) = 
J (i^n$(r, n)i^J^(n), and ai{r) was defined previously (Eq.[5]). We can instead use isocurvature 
fluctuations or a mixture of the two. Equation (Al) suggests that a;(r) is a transfer function 
projecting $(x) onto Ar(n) through the integral over the line of sight. Since a;(r) is just a math- 
ematical function, we pre-compute and store it for a given cosmology, reducing the computational 
time of a batch of simulations. We can thus use or extend equation (Al) to compute Ar(n) for 
generic primordial fluctuations. 

We simulate CMB sky maps using a non-Gaussian model of the form in equation (1) as follows. 
(1) We generate 'I>L(k) as a Gaussian random field in Fourier space for a given initial power spectrum 
P{k), and transform it back to real space to obtain <I>l(x). (2) We transform from Cartesian to 
spherical coordinates to obtain ^i^{r,n), compute its harmonic coefficients ^/^(r), and obtain a 
temperature map of the Gaussian part Ar$(n) by integrating equation (Al). (3) We repeat this 
procedure for ~ ^x'^ I c?^x$l(x) to obtain a temperature map of the non-Gaussian part 
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Ar$2(n). (4) By combining these two temperature maps, we obtain non-Gaussian sky maps for 
any values of /nl, 

AT(n) = AT$(n) + /NLAr$2 (n). (A2) 

We do not need to run many simulations individually for different values of /nl, but run only twice 
to obtain AT$(n) and Ar,j,2(n) for a given initial random number seed. Also, we can combine 
Ar$(n) for one seed with Ar$2(n) for the other to make realizations for a particular kind of 
two-field inflation models. We can apply the same procedure to isocurvature fluctuations with or 
without $(x) correlations. 

We need the simulation box of the size of the present-day cosmic horizon size L^ox = 2cro, 
where tq is the present-day conformal time. For example, Lbox ~ 20 Gpc is needed for a flat 
universe with = 0.3, whereas we need spatial resolution of at least ~ 20 Mpc to resolve 
the last-scattering surface accurately. Prom this constraint the number of grid points is at least 
-^grid = 1024^, and the required amount of physical memory to store $(x) is at least 4.3 GB. 
Moreover, when we simulate a sky map having 786 432 pixels at nside = 256, we need 1.6 GB to 
store a field in spherical coordinate $(r, n), where the number of r evaluated for Ngs-id = 1024^ 
is 512. Since our algorithm for transforming Cartesian into spherical coordinates requires another 
1.6 GB, in total we need at least 7.5 GB of physical memory to simulate one sky map. 

We have generated 300 realizations of non-Gaussian sky maps with A^grid = 1024^ and nside = 
256. It takes 3 hours on 1 processor of SGI Origin 300 to simulate Ar$(n) and Ar$2(n). We 
have used 6 processors to simulate 300 maps in one week. Figure 7 shows the one-point probability 
density function (PDF) of temperature fluctuations measured from simulated non-Gaussian maps 
(without noise and beam smearing) compared with the r.m.s. scatter of Gaussian realizations. We 
find it difficult for the PDF alone to distinguish non-Gaussian maps of |/nl| < 500 from Gaussian 
maps, whereas the cubic statistic -Sprim (eq-[8]) can easily detect /nl = 100 in the same data sets. 

We measure /nl on the simulated maps using S'prim to see if it can accurately recover /nl- 
Similar tests show the Minkowski functionals to be unbiased and able to discriminate different /nl 
values at levels consistent with the quoted uncertainties. We also measure the point-source angular 
bispectrum ftgrc to see if it returns null values as the simulations do not contain point sources. We 
have included noise properties and window functions in the simulations. Figure 8 shows histograms 
of /nl and 6src measured from 300 simulated maps of /nl = 100 (solid lines) and /nl = (dashed 
lines). Our statistics find correct values for /nl and find null values for ftgrc; thus, our statistics 
are unbiased, and /nl and bg^c 

are orthogonal to each other as pointed out by Komatsu & Spergel 

(2001). 

B. POWER OF THE POINT-SOURCE BISPECTRUM 

In this appendix, we test our estimator for fegrc and Cgrc using simulated Q-band maps of point 
sources, CMB, and detector noise. The 44 GHz source count model of T98 was used to generate 
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Fig. 7. — One-point PDF of temperature fluctuations measured from simulated non-Gaussian maps 
(noise and beam smearing are not included). From the top-left to the bottom-right panel the solid 
lines show the PDF for /nl = 100, 500, 1000, and 3000, while the dashed lines enclose the r.m.s. 
scatter of Gaussian realizations (i.e., /nl = 0). 
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Fig. 8. — The distribution of the non-hnear coupling parameter /nl (tlic left panel) and the point- 
source bispectrum bsrc (the right panel) measured from 300 simulated realizations of non-Gaussian 
maps for /nl = 100 (solid line) and /nl = (dashed line). The simulations include noise properties 
and window functions of the WMAP 1-year data, but do not include point sources. 
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the source populations. The total source count in each realization was fixed to 9043, the number 
predicted by T98 to lie between S'min = O.lJy and Smax = 10 Jy. By generating uniform deviates 
■u G (0, 1) and transforming to flux S via 

^ iV(> 5n,in) - Ni> S) 

N{> S^in) - N{> S^^) ' 

we obtain the desired spectrum. The sources were distributed evenly over the sky and convolved 
with a Gaussian profile approximating the Q-band beam. Flux was converted to peak brightness 
using the values in Table 8 of Page et al. (2003b). The CMB and noise realizations were not varied 
between realizations. The goal in this appendix is to prove that our estimator for 6src works well 
and is very powerful in estimating Cgrc- 

The left panel of Figure 9 compares the measured 6src from simulated maps with the expec- 
tations of the simulations. Black, dark-gray, and light-gray indicate three different realizations of 
point sources. The measurements agree well with the expectations at < 1.75 Jy. They however 
show significant scatter at Sc > 1.75 Jy, because our filter for computing bgic [eq. (21)] does not 

include contribution from Cgrc to Ci, making the filter less optimal in the limit of "too many" un- 
masked point sources. We can see from the figure that Cgrc at > 2 Jy is comparable to or larger 
than the noise power spectrum for Q band, 54 x 10"'^ //K^ sr. 

Fortunately this is not a problem in practice, as we can detect and mask those bright sources 
which contribute significantly to Q. The residual point sources that we cannot detect (therefore 
we want to quantify using ftgrc) should be hidden in the noise having only a small contribution to 
Q. In this faint-source regime 6src works well in measuring the amplitude of residual point sources, 
offering a promising way for estimating Cgrc- The right panel of Figure 9 compares Cgrc estimated 
from bsTc [Eq. (17)] with the expectations. The agreement is good for Sc < 1.75 Jy, proving that 
estimates of Cgrc from 6src are unbiased and powerful. Since 6src measures Cgrc directly, we can use 
it for any CMB experiments which suffer from the effect of residual point sources. While we have 
considered the bispectrum only here, the forth order moment may also be used to increase our 
sensitivity to the point-source non-Gaussianity (Pierpaoli 2003). 



(Bl) 



C. THE ANGULAR BISPECTRUM FROM A POTENTIAL STEP 

A scalar-field potential V{(f)) with features can generate large non-Gaussian fluctuations in 
CMB by breaking the slow-roll conditions at the location of the features (Kofman et al. 1991; 
Wang &; Kamionkowski 2000). We estimate the impact of the features by using a scale-dependent 
/nl, 

which is calculated from a non-linear transformation between the curvature perturbations in the 
comoving gauge and the scalar-fleld fluctuations in the spatially flat gauge (Salopek & Bond 1990, 
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Fig. 9. — Testing the estimator for the reduced point-source bispectrum fegrc [Eq.(22)]. The left 
panel shows bgrc measured from a simulated map including point sources and properties of the 
WMAP sky map at Q band, as a function of flux cut Sc (filled circles). Black, dark-gray, and 
light-gray indicate three different realizations of point sources. The solid line is the expectation 
from the input source number counts in the point-source simulation. The right panel compares the 
power spectrum Cgrc estimated from 6src with the expectation. The error bars are not independent, 
because the distribution is cumulative. The behavior for > 2 Jy shows the cumulative effect of 
sources with brightness comparable to the instrument noise (see text in Appendix B). 



-30- 



1991). This expression does not assume the slow-roll conditions. Although this expression does not 
include all effects contributing to /nl during inflation driven by a single field (Maldacena 2002), 
we assume that an order-of-magnitude estimate can still be obtained. 

A sharp feature in V{(f)) at (pf produces a significantly scale-dependent /nl(i;^) near (pf through 
the derivatives of H in equation (CI). We illustrate the effects of the steps using the potential 
features proposed by Adams et al. (2001), 

1 " 



1 + ctanh 



d 



(C2) 



which has a step in V^cf)) at (/){ with the height c and the slope dr^. Adams ct al. (1997) show 
that the steps are created by a class of supergravity models in which symmetry-breaking phase 
transitions of many fields in flat directions gravitationally coupled to continuously generate steps 
in V{<p) every 10 — 15 e-folds, giving a chance for a step to exist within the observable region of 
V{cf>). 

It is instructive to evaluate equation (CI) combined with equation (C2) in the slow-roll limit, 
In H/d(f)'^ ~ ^d^lnV/dcl)^. For |c| < 1, one obtains 

where x = (/){)/ d. The first term corresponds to a standard, nearly scale-independent prediction 
giving 7.4 x 10^'^ at cj) = Smpianck) while the second term reveals a significant scale-dependence. 
The function tanhx/ cosh^ a; is a symmetric odd function about x = with extrema of ±0.385 at 
X ~ ±0.66. The picture is the following: as (p rolls down V{(p) from a positive x > 0.66, (p gets 
accelerated at x ~ 0.66, reaches constant velocity at a; = 0, decelerates at x ~ —0.66, and finally 
reaches slow roll at a; < —0.66. The ratio of the second term in equation (C3) to the first at the 
extrema is ±0.385c(^/(i)^. For example, c = 0.02 and (p/d = 300 (i.e., d = O.Olmpianck) make the 
amplitude of the second term 700 times larger than the first, giving |/nl| — 5 at the extrema. 
Despite the slow-roll conditions having a tendency to underestimate /nl, it is possible to obtain 
|/nl| > 1- Neglecting the first term in equation (C3) and converting (p for k, one obtains 

. 5c , ,,, 5c tanhxfc 

" 24;^^«*^P(^) " 24vrGd^cosh^x. ' ^^'^ 

where x/- — d~^ {d(p/d In k){{k/kf — 1) = d^^{(p/ H)i{k/ki — 1) for k — kf <^ fcf. The slow-roll 
approximation gives Xk — {'^nGcpid)^^{k/ki — 1). Finally, following the method of Komatsu & 
Spergel (2001), we obtain the reduced bispectrum of a potential step model, bfj^i^, as 

= 2 (24^) /"^'^^ [PiAm,{r)ot^:^{r) + (2 permutations)] , (C5) 
where /3z(r) is given by equation (6), and 

af P(r) ='^j k^dkh,tepik)gTiik)jiikr). (C6) 
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The amplitude is thus proportional to c/d^: a bigger (larger c) and steeper (smaller d) step gives 
a larger bispcctrum. The steepness affects the amplitude more, because the non-Gaussianity is 
generated by breaking the slow-roll conditions. 

Since bf^'i^i-^ linearly scales as c for a fixed d, we can fit for c by using exactly the same method as 
for the scale-independent /nl, but with ai{r) in equation (3) replaced by af''^'^{r). The exact form 
of the fitting parameter in the slow-roll limit is 5c/(247rGci^). A reason for the similarity between 
the two models in methods for the measurement is explained as follows. Komatsu et al. (2003) have 
shown that B{h,r) [Eq.(4)] is a Wiener-filtered, reconstructed map of the primordial fluctuations 
$(n, r). Our cubic statistic [Eq.(7)] effectively measures the skewness of the reconstructed $ field, 
maximizing the sensitivity to the primordial non-Gaussianity. One of the three maps comprising 
the cubic statistic is however not B{h, r), but A{h, r) given by equation (3). This map defines what 
kind of non-Gaussianity to look for, or more detailed form of the bispectrum. For the potential 
step case, Astep(n, r) made of af^^{r) picks up the location of the step to measure 5c/(247rG(i^) 
near k{, while for the form in equation (1), A{n,r) explores all scales on equal footing to measure 
the scale-independent /nl- 

The distinct features in k space are often smeared out in I space via the projection. This effect is 
estimated from equation (C6) as follows. The function /istep(fc) near kf is accurately approximated 
by hstcp{k) — 0.385 sin(2xfc), which has a period of Ak = 47r^G0fdA;f. On the other hand, the 
radiation transfer function gTi{k) behaves as ji{kr^) where is the comoving distance to the photon 
decoupling epoch, and gTi{k)ji{kr) behaves as jf{kr^) (the integral is very small when r 7^ r*). The 
oscillation period of this part is thus AA; = 7r/r* for /cr* > Z. A ratio of the period of /istep(fc) to 
that of gTi{k)ji{kr) is then estimated as AiTG(f)idr^ki ~ (Zf/3)(d/0.01mpianck)(<?^f/3mpianck)) where 
Zf = fcfr* is the angular wave number of the location of the step. We thus find that hstep{k) oscillates 
much more slowly than the rest of the integrand in equation (C6) for Zf » 1. 

What does it mean? It means that the results would look as if there were two distinct regions 
in I space where /nl is very large: a positive /nl is found at I < l{ and a negative one at Z > Zf. 
The estimated location is l/l{ ~ 1 ± 0.66(47rG(/)fd) ~ 1 it 0.2(d/0.01mpianck)(</'f/3?Ti-pianck); thus, 
the positive and negative regions are separated in I by only 40%, making the detection difficult 
when many Z modes are combined to improve the signal-to-noise ratio. The two extrema would 
cancel out to give only small signals. In other words, it is still possible that non-Gaussianity from 
a potential step is "hidden" in our measurements shown in Figure 1. Note that the cancellation 
occurs because of the point symmetry of Zistep(fc) about k = kf. If the function has a knee instead 
of a step, then the cancellation does not occur and there would be a single region in I space where 
|/nl| is large (Wang & Kamionkowski 2000). Note that our estimate in this Appendix was based 
upon equation (C3), which uses the slow-roll approximations. While instructive, since the slow-roll 
approximations break down near the features, our estimate may not be very accrate. One needs 
to integrate the equation of motion of the scalar field to evaluate equation (CI) for more accurate 
estimations of the effect. 
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